#include "random.hh"

random_index_t::random_index_t(int K): kmax(K)
{
#ifdef DEBUG
  assert( K > 1 );
#endif
  RandGen = new boost::random::mt19937(time(0));
}

random_index_t::~random_index_t()
{
  delete RandGen;
}

int
random_index_t::operator () (const sp_stat_vec_t& pivec)
  const
{
  boost::random::uniform_int_distribution<> dist(1, 100000);
  double tot = sum( pivec, 1 );
  double u = dist(*RandGen) / (100000.) * tot;
  double cumsum = 0.;
  for(sp_stat_vec_t::iterator_t ki=pivec.begin_nonzero();
      ki!=pivec.end_nonzero(); ++ki)
    {
      int k = ki.dim();
      double p = ki.value();
      cumsum += p;
      if( u <= cumsum )
        return k;
    }
  return (kmax-1);
}

int
random_index_t::operator () ()
  const
{
  boost::random::uniform_int_distribution<> dist(0,kmax-1);
  return dist(*RandGen);
}

int
random_index_t::operator () (int K)
  const
{
#ifdef DEBUG
  assert( K > 1 );
#endif
  boost::random::uniform_int_distribution<> dist(0,K-1);
  return dist(*RandGen);
}
